library(tidyverse)
theme_set(theme_classic())

Le package ggmap

ggmap est un package permettant de faire de la cartographie. On pourra consulter cet article pour plus de détails.

ggmap permet de récupérer facilement des fonds de carte. Par exemple :

library(ggmap)
us <- c(left = -125, bottom = 25.75, right = -67, top = 49)
map <- get_stamenmap(us, zoom = 5, maptype = "toner-lite")
ggmap(map)

On peut utiliser dplyr pour les opérations

get_stamenmap(us, zoom = 5, maptype = "toner-lite") %>% ggmap()

Pour l’Europe on fait

europe <- c(left = -12, bottom = 35, right = 30, top = 63)
get_stamenmap(europe, zoom = 5,"toner-lite") %>% ggmap()

On peut également changer le fond de carte

get_stamenmap(europe, zoom = 5,"toner-background") %>% ggmap()

Pour la france, on aura

fr <- c(left = -6, bottom = 41, right = 10, top = 52)
get_stamenmap(fr, zoom = 5,"toner-lite") %>% ggmap()

La fonction geocode qui permettait de récupérer des latitudes et longitudes n’est plus gratuite. On propose d’utiliser la fonction suivante :

if (!(require(jsonlite))) install.packages("jsonlite")
geocodeGratuit <- function(adresses){
# adresses est un vecteur contenant toutes les adresses sous forme de chaine de caracteres
  nominatim_osm <- function(address = NULL){
    ## details: http://wiki.openstreetmap.org/wiki/Nominatim
    ## fonction nominatim_osm proposée par D.Kisler
    if(suppressWarnings(is.null(address)))  return(data.frame())
    tryCatch(
      d <- jsonlite::fromJSON(
        gsub('\\@addr\\@', gsub('\\s+', '\\%20', address),
             'http://nominatim.openstreetmap.org/search/@addr@?format=json&addressdetails=0&limit=1')
      ), error = function(c) return(data.frame())
    )
    if(length(d) == 0) return(data.frame())
    return(c(as.numeric(d$lon), as.numeric(d$lat)))
  }
  tableau <- t(sapply(adresses,nominatim_osm))
  colnames(tableau) <- c("lon","lat")
  return(tableau)
}

Cette fonction permet de récupérer les latitudes et longitudes de lieux à spécifier :

geocodeGratuit("the white house")

##                       lon     lat
## the white house -77.03655 38.8977

geocodeGratuit("Paris")

##            lon     lat
## Paris 2.351462 48.8567

geocodeGratuit("Rennes")

##             lon      lat
## Rennes -1.68002 48.11134

Exercice 1

Récupérer les latitudes et longitudes de Paris, Lyon et Marseille et représenter ces 3 villes sur une carte de la france.

V <- c("Paris","Lyon","Marseille")
A <- geocodeGratuit(V)
A <- A %>% as_tibble() %>% mutate(Villes=V)
fr <- c(left = -6, bottom = 41, right = 10, top = 52)
fond <- get_stamenmap(fr, zoom = 5,"toner-lite") 
ggmap(fond)+geom_point(data=A,aes(x=lon,y=lat),color="red")

Le fichier villes_fr.csv contient les populations des 30 plus grandes villes de france. Représenter à l’aide d’un point les 30 plus grandes villes de france. On fera varier la taille du point en fonction de la population en 2014.

df <- read_csv("villes_fr.csv")
df$Commune <- as.character(df$Commune)
df$Commune[10] <- "Lille"
#coord <- geocodeGratuit(as.character(df$Commune)) %>% as_tibble()
coord <- read.csv("coord_ville_exo1.csv")
df1 <- bind_cols(df,coord)
df1 <- cbind.data.frame(df,coord)
ggmap(fond)+geom_point(data=df1,aes(x=lon,y=lat),color="red")

ggmap(fond)+geom_point(data=df1,aes(x=lon,y=lat,size=`2014`),color="red")

Enfin, il est bien entendu possible de faire des cartes en ayant à disposition des coordonnées de frontières. Ces coordonnées peuvent se trouver sur le net ou dans certains packages R tels que maps :

library(maps)
france <- map(database="france")

names(france)

## [1] "x"     "y"     "range" "names"

On peut alors utiliser ggplot :

fr1 <- tibble(lon=france$x,lat=france$y)
ggplot(fr1)+geom_path(aes(y=lat,x=lon),size=1)+theme_classic()

On pourra supprimer les axes avec

ggplot(fr1)+geom_path(aes(y=lat,x=lon),size=1)+theme_void()

Le package sf

L’approche précédente pour représenter les département en ggplot n’est pas vraiment (ou pas du tout) satisfaisante. En effet, nous représentons les latitudes et longitudes sur un repère cartésien 2D. On rappelle que pour visualiser une carte en 2D il faut réaliser des projections. Ces projections peuvent bien entendu se calculer la main mais les packages spécifiques (comme ggmap) vont les faire automatiquement. Nous proposons ici de présenter brièvement le package sf qui va nous permettre de créer des cartes “avancées”. On pourra trouver de la documentation sur ce package aux url suivantes :

https://statnmap.com/fr/2018-07-14-initiation-a-la-cartographie-avec-sf-et-compagnie/
dans les vignettes sur la page du cran de ce package : https://cran.r-project.org/web/packages/sf/index.html

Ce package propose de définir un nouveau format sf adapté à la cartographie. Regardons par exemple l’objet nc

library(sf)
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc

Cet objet contient des informations sur les morts subites de nourissons dans des villes de Caroline du Nord. On remarque la présence d’une colonne geometry qui permettra de délimiter ces villes à l’aide de polygones. L’objet défini, il est facile de visualiser la carte avec un plot classique

plot(st_geometry(nc))

ou avec ggplot

ggplot(nc)+geom_sf()

On peut également ajouter des labels dans les polygones

ggplot(nc[1:3,]) +
   geom_sf(aes(fill = AREA)) + 
   geom_sf_label(aes(label = NAME))

On peut représenter les villes de la table avec

#coord.ville <- matrix(0,ncol=2,nrow=length(nc$NAME))
#for (i in 1:length(nc$NAME)){
#  i
#  coord.ville[i,] <- geocodeGratuit(as.character(nc$NAME[i]))
#}
coord.ville <- read.csv("coord.ville_sf.csv")
coord.ville <- as.data.frame(coord.ville)
names(coord.ville) <- c("lon","lat")

#coord.ville <- geocodeGratuit(nc$NAME)
coord.ville1 <- coord.ville %>%  
  filter(lon<=-77 & lon>=-85 & lat>=33 & lat<=37) %>% 
  as.matrix() %>% st_multipoint() %>% st_geometry()
st_crs(coord.ville1) <- 4326 
ggplot(nc)+geom_sf()+geom_sf(data=coord.ville1)

Il est souvent préférable de considérer un objet sf sous la forme d’un dataframe ou tibble. On peut par exemple le faire pour l’objet coord.ville1 avec

coord.ville2 <- st_cast(coord.ville1, to = "POINT")
nom_ville <- data.frame(coord.ville[,1:2]) %>% 
  mutate(NAME=nc$NAME) %>% 
  filter(lon<=-77 & lon>=-85 & lat>=33 & lat<=37)
df <- data.frame(Ville=nom_ville)
st_geometry(df) <- coord.ville2
st_crs(df) <- 4326
ggplot(nc)+geom_sf()+geom_sf(data=df)

Exercice 2

Nous nous servons de la carte GEOFLAR proposée par l’Institut Géographique National pour récupérer un fond de carte contenant les frontières des départements français. Cette carte est disponible sur le site http: //professionnels.ign.fr/ au format shapefile, elle se trouve dans l’archive dpt.zip. Il faut décompresser pour reproduire la carte. Grâce au package sf, cette carte, contenue dans la série de fichiers departement du répertoire dpt, peut être importée dans un objet R :

dpt <- read_sf("dpt")
ggplot(dpt) + geom_sf()

Refaire la carte de l’exercice 1 sur ce fond de carte.

df <- read_csv("villes_fr.csv")
df$Commune <- as.character(df$Commune)
df$Commune[10] <- "Lille"
coord <- geocodeGratuit(as.character(df$Commune)) %>% as_tibble()
df1 <- bind_cols(df,coord)

coord.ville1 <- data.frame(df1[,14:15]) %>% 
  as.matrix() %>% st_multipoint() %>% st_geometry()

coord.ville2 <- st_cast(coord.ville1, to = "POINT")

coord.ville1

## Geometry set for 1 feature 
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: -4.486009 ymin: 42.69853 xmax: 7.750713 ymax: 50.63657
## epsg (SRID):    NA
## proj4string:    NA

## MULTIPOINT (2.351462 48.8567, 5.369953 43.29617...

coord.ville2

## Geometry set for 30 features 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -4.486009 ymin: 42.69853 xmax: 7.750713 ymax: 50.63657
## epsg (SRID):    NA
## proj4string:    NA
## First 5 geometries:

## POINT (2.351462 48.8567)

## POINT (5.369953 43.29617)

## POINT (4.832011 45.75781)

## POINT (1.444247 43.60446)

## POINT (7.268391 43.70094)

st_geometry(df1) <- coord.ville2
st_crs(df1) <- 4326
ggplot(dpt)+geom_sf(fill="white")+
  geom_sf(data=df1,aes(size=`2014`),color="red")+theme_void()

Exercice 3

Nous souhaitons visualiser graphiquement les différences de taux de chômage par département entre deux années. Pour cela, nous disposons de chaque taux mesuré aux premiers trimestres des années 2006 et 2011 (variables TCHOMB1T06, TCHOMB1T11) qui se trouvent dans le jeu de données tauxchomage.csv

Importer le jeu de données.

chomage <- read_delim("tauxchomage.csv",delim=";")

Faire la jointure de cette table avec celle des départements. On pourra utiliser inner_join.

dpt <- read_sf("dpt")
dpt2 <- inner_join(dpt,chomage,by="CODE_DEPT")

Comparer les taux de chomage en 2006 et 2011.

dpt3 <- dpt2 %>% select(A2006=TCHOMB1T06,A2011=TCHOMB1T11,geometry) %>%
  gather("Annee","TxChomage",-geometry)

ggplot(dpt3) + aes(fill = TxChomage)+geom_sf() +
  facet_wrap(~Annee, nrow = 1)+
  scale_fill_gradient(low="white",high="brown")+theme_bw()

Challenge 1 : carte des températures avec ggmap

On importe d’abord les données contenant des températuers observées ainsi que des coordonées de stations.

Les fichiers synop.2018011215.csv et postesSynop.csv contiennent les températures (colonne t) du 12 janvier 2018 à 15h et les localisations de 59 stations françaises. Concaténer les deux tables en utilisant comme clef le numéro de la station. On pourra lire les tables avec read_delim et faire la jointure avec inner_join.

donnees <- read_delim("synop.2018011215.csv",delim=";")
station <- read_delim("postesSynop.csv",delim=";")
temp <- donnees %>% select(numer_sta,t)
names(temp)[1] <- c("ID")
D <- inner_join(temp, station, by = c("ID"))

Faire de même avec les données actuelles. On pourra requêter directement sur le site des données publiques de meteo france en utilisant (par exemple) la fonction fread du package data.table ou la fonction read_delim du package readr. Convertir les degrés Kelvin en degrés Celsius.

donnees <- read_delim("https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/synop.2020020115.csv",delim=";")
station <- read_delim("https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv",delim=";")
donnees$t <- donnees$t-273.15 #on passe en degrés celcius
temp <- donnees %>% select(numer_sta,t)
names(temp)[1] <- c("ID")
D <- inner_join(temp, station, by = c("ID"))

Eliminer les station d’outre mer (on pourra conserver uniquement les stations qui ont une longitude entre -20 et 25). On appellera ce tableau station1. Visualiser les stations sur la carte contenant les frontières des départements français.

station1 <- D %>% filter(Longitude<25 & Longitude>-20)
station4326 <- st_multipoint(as.matrix(station1[,5:4])) %>% st_geometry()
st_crs(station4326) <- 4326
ggplot(dpt) + geom_sf()+geom_sf(data=station4326)

Créer un dataframe au format sf qui contient les températures des stations ainsi que leurs coordonnées dans la colonne geometry. On pourra commencer avec

station2 <- station1 %>% select(Longitude,Latitude) %>% as.matrix() %>% st_multipoint() %>% st_geometry()
st_crs(station2) <- 4326
station2 <- st_cast(station2, to = "POINT")

df <- data.frame(temp=station1$t)
st_geometry(df) <- station2

Représenter les stations sur une carte de france. On pourra mettre un point de couleur différente en fonction de la température.

ggplot(dpt) + geom_sf(fill="white")+
  geom_sf(data=df,aes(color=temp),size=2)+
  scale_color_continuous(low="yellow",high="red")

On obtient les cordonnées des centroïdes des départements à l’aide de

centro <- st_centroid(dpt$geometry) 
centro <- st_transform(centro,crs=4326)

On déduit les distances entre ces centroïdes et les stations avec (df étant la table sf obtenue à la question 3).

DD <- st_distance(df,centro)

Prédire la température de chaque département à l’aide de la règle du 1 plus proche voisin.

NN <- apply(DD,2,order)[1,]
t_prev <- station1[NN,2]

Colorier les départements en fonction de la température prédite dans le département. On pourra faire varier le dégradé de couleur du jaune (pour les faibles températures) au rouge (pour les fortes).

dpt1 <- dpt %>% mutate(t_prev=as.matrix(t_prev))
ggplot(dpt1) + geom_sf(aes(fill=t_prev)) +
  scale_fill_continuous(low="yellow",high="red")

Le package leaflet

Leaflet est un package permettant de faire de la cartographie intéractive. On pourra consulter un descriptif synthétique ici. On peut par exemple obtenir un fond de carte de la terre avec

library(leaflet)
m <- leaflet()
m <- addTiles(m)
m

Ou de manière équivalente avec dplyr

leaflet() %>% addTiles()

On dispose de plusieurs styles de fonds de cartes (quelques exemples ici) :

Paris <- geocodeGratuit("paris")
#m2 <- leaflet() %>% setView(lng = Paris[1], lat = Paris[2], zoom = 12) %>% 
#  addTiles()
#m2 %>% addProviderTiles("Stamen.Toner")
#m2 %>% addProviderTiles("OpenStreetMap.BlackAndWhite")
#m2 %>% addProviderTiles("Thunderforest.Transport")

On pourra placer des markers sur la carte avec des commandes telles que addMarkers ou addCircles…

data(quakes)
leaflet(data = quakes[1:20,]) %>% addTiles() %>%
  addMarkers(~long, ~lat, popup = ~as.character(mag))

Les popups permettent d’ajouter des informations auxiliaires lorsqu’on clique sur un marker dans la carte :

content <- paste(sep = "<br/>",
  "<b><a href='http://www.samurainoodle.com'>Samurai Noodle</a></b>",
  "606 5th Ave. S",
  "Seattle, WA 98138"
)

leaflet() %>% addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
  addPopups(-122.327298, 47.597131, content,
    options = popupOptions(closeButton = FALSE)
  )

Exercice 4

Placer un popup représentant l’Université Rennes 2 (Campus Villejean). On ajoutera un lien renvoyant sur le sité de l’Université.

R2 <- geocodeGratuit("Universite Rennes 2 Villejean") %>% as_tibble()
info <- paste(sep = "<br/>",
  "<b><a href='https://www.univ-rennes2.fr'>Universite Rennes 2</a></b>",
  "Campus Villejean")


leaflet() %>% addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%  
  addPopups(R2[1]$lon, R2[2]$lat, info,options = popupOptions(closeButton = FALSE))

Challenge 2 : Velib Paris

Plusieurs villes dans le monde ont accepté de mettre en ligne les données sur l’occupation des stations velib dans leur ville. Ces données sont facilement accessibles et mises à jour en temps réel. On dispose généralement de la taille et la localisation des stations, la proportion de vélos disponibles… Il est possible de requêter (entre autres) :

sur les données Decaux
sur Open data Paris
sur vlstats pour des données mensuelles ou historiques ou encore sur Velib pour obtenir des fichiers qui sont rafraichis régulièrement.

Récupérer les données actuelles de velib disponibles pour la ville de Paris : https://opendata.paris.fr/explore/dataset/velib-disponibilite-en-temps-reel/information/. On pourra utiliser la fonction url dans read.csv.

sta.Paris <- read.csv(url("https://opendata.paris.fr/explore/dataset/velib-disponibilite-en-temps-reel/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true"),sep=";")

Décrire les variables du jeu de données.

Nous avons de l’information sur la disponibilité, le remplissage… de stations velib parisiennes.

Créer les une variable latitude et une variable longitude à partir de la variable geo.

sta.Paris1 <- sta.Paris %>% separate(geo,into=c("lat","lon"),sep=",") %>% 
  mutate(lat=as.numeric(lat),lon=as.numeric(lon))

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [354,
## 717].

Visualiser les positions des stations sur une carte leaflet.

map.velib1 <- leaflet(data = sta.Paris1) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
  addCircleMarkers(~ lon, ~ lat,radius=3,
    stroke = FALSE, fillOpacity = 0.5,color="red"
  )

map.velib1

Ajouter un popup qui permet de connaitre le nombre de vélos disponibles (électriques+mécanique) quand on clique sur la station (on pourra utiliser l’option popup dans la fonction addCircleMarkers).

sta.Paris2 <- sta.Paris1 %>% 
  mutate(nb_velo_dispo=Nombre.de.vélo.mécanique+Nombre.vélo.électrique)
map.velib2 <- leaflet(data = sta.Paris2) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
  addCircleMarkers(~ lon, ~ lat,radius=3,stroke = FALSE, 
                   fillOpacity = 0.7,color="red", 
                   popup = ~ sprintf("<b> Vélos dispos: %s</b>",as.character(nb_velo_dispo)))

#ou sans sprintf

map.velib2 <- leaflet(data = sta.Paris2) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>% 
  addCircleMarkers(~ lon, ~ lat,radius=3,stroke = FALSE, fillOpacity = 0.7,color="red", 
                   popup = ~ paste("Vélos dispos :",as.character(nb_velo_dispo)))

map.velib2

Ajouter la nom de la station dans le popup.

map.velib3 <- leaflet(data = sta.Paris2) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
  addCircleMarkers(~ lon, ~ lat,radius=3,stroke = FALSE, fillOpacity = 0.7,color="red", 
                   popup = ~ paste(as.character(Nom.de.la.station),", Vélos dispos :",as.character(nb_velo_dispo),sep=""))

## Warning in validateCoords(lng, lat, funcName): Data contains 2 rows with
## either missing or invalid lat/lon values and will be ignored

map.velib3

Faire de même en utilisant des couleurs différentes en fonction de la proportion de vélos disponibles dans la station. On pourra utiliser les palettes de couleur

ColorPal1 <- colorNumeric(scales::seq_gradient_pal(low = "#132B43", high = "#56B1F7", space = "Lab"), domain = c(0,1))
ColorPal2 <- colorNumeric(scales::seq_gradient_pal(low = "red", high = "black", space = "Lab"), domain = c(0,1))

map.velib4 <- leaflet(data = sta.Paris2) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
  addCircleMarkers(~ lon, ~ lat,radius=3,stroke = FALSE, fillOpacity = 0.7,
                   color=~ColorPal1(nb_velo_dispo/Nombres.de.bornes.en.station), 
                   popup = ~ paste(as.character(Nom.de.la.station),", Vélos dispos :",as.character(nb_velo_dispo),sep=""))

map.velib4

map.velib5 <- leaflet(data = sta.Paris2) %>% 
  addTiles(urlTemplate = "http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
  addCircleMarkers(~ lon, ~ lat,radius=3,stroke = FALSE, fillOpacity = 0.7,
                   color=~ColorPal2(nb_velo_dispo/Nombres.de.bornes.en.station), 
                   popup = ~ paste(as.character(Nom.de.la.station),", Vélos dispos :",as.character(nb_velo_dispo),sep=""))

map.velib5

Exercice 5

Refaire la carte des températures du premier challenge en utilisant leaflet. On utilisera la table construite dans le challenge 1 et la fonction addPolygons. On pourra également ajouter un popup qui permet de visualiser le nom du département ainsi que la température prévue lorsqu’on clique dessus.

factpal <- colorFactor(topo.colors(5), dpt1$t_prev)
pal <- colorNumeric(palette = c("inferno"),domain = dpt1$t_prev)
dpt2 <- st_transform(dpt1, crs = 4326)
dpt2$t_prev <- round(dpt2$t_prev)
m <- leaflet() %>% addTiles() %>% 
  addPolygons(data = dpt2,color=~pal(t_prev),fillOpacity = 0.6, 
              stroke = TRUE,weight=1,
              popup=~paste(as.character(NOM_DEPT),as.character(t_prev),sep=" : ")) %>% 
  addLayersControl(options=layersControlOptions(collapsed = FALSE))

## Warning in pal(t_prev): Some values were outside the color scale and will
## be treated as NA

## Warning in pal(t_prev): Some values were outside the color scale and will
## be treated as NA

Cartographie : ggmap, sf et leaflet